Constrained-path quantum Monte Carlo simulations of the zero-temperature, 
disordered two-dimensional Hubbard model 
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We study the effects of disorder on long-range antiferro- 
magnetic correlations in the half-filled, two dimensional, re- 
pulsive Hubbard model at T = 0. A mean field approach is 
first employed to gain a qualitative picture of the physics and 
to guide our choice for a trial wave function in a constrained 
path quantum Monte Carlo (CPQMC) method that allows for 
a more accurate treatment of correlations. Within the mean 
field calculation, we observe both Anderson and Mott insulat- 
ing antiferromagnetic phases. There are transitions to a para- 
magnet only for relatively weak coupling, U < 2t in the case 
of bond disorder, and U < 4t in the case of on-site disorder. 
Using ground state CPQMC we demonstrate that this mean 
field approach significantly overestimates magnetic order. For 
U = 4t, we find a critical bond disorder of 14 ~ (1.6 ± 0.4)t 
even though within mean field theory no paramagnetic phase 
is found for this value of the interaction. In the site disordered 
case, we find a critical disorder of Vc ~ (5.0±0.5)f at U — At. 
PACS numbers: 74.20.-z, 74.20.Mn, 74.25.Dw 



I. INTRODUCTION 

The Hubbard Hamiltonian encapsulates many of the 
most interesting qualitative many-body effects in corre- 
lated fermion systems, notably the possibility of ordering 
of electron spins and the appearance of insulating states 
in systems with partially filled electronic bands. While 
the original model is translationally invariant, the intro- 
duction of disorder can alter the magnetic and charge 
correlations in fundamental ways. [QH] Experimental sys- 
tems whose qualitative physics appears to involve the in- 
terplay of interactions and randomness, possibly modeled 
by the Hubbard Hamiltonian, include doped semiconduc- 
tors H,^, thin superconducting films |Q,D, and silicon 
mctal-oxide-semiconductor field-effect transistors. Q 

The determinant quantum Monte Carlo (DQMC) 
method has been a useful tool to simulate the Hubbard 
Hamiltonian, ||^,^ but its application has been limited by 
the impossibility of reaching low temperatures in many 
cases of interest, a problem which arises when the "Boltz- 
mann weight" for the fermion system becomes negative. 

An approach for dealing with this "sign problem" in 
the ground state is the constrained path quantum Monte 
Carlo (CPQMC) method. @ In CPQMC, the sign prob- 
lem is treated by imposing, in the space of Slater deter- 
minants, a boundary condition based on an input trial 



wave function. 

The CPQMC approach has not previously been used 
in models with quenched randomness. In this paper, we 
apply CPQMC to the disordered, two-dimensional (2D), 
"Anderson-Hubbard" Hamiltonian, 
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Here c\c (c;^) are operators that destroy (create) elec- 
trons of spin a on site i of a 2D square lattice of size 
LF' = N.U is the on-site repulsion, /i and ti are the chem- 
ical potential and random site energies, respectively, and 
tij is the (random) hopping energy. Random on-site ener- 
gies are chosen uniformly on [—Vs/2,-\-Vs/2], and are 
chosen uniformly on [t — Vt/2,t-\- Vt /2] , where Vs and Vt 
are parameters that set the disorder strength. We will 
choose t — 1 to set the scale of energy, and focus our at- 
tention on the case when the lattice is half-filled {n) — 1. 

In the absence of disorder, the half-filled Hubbard 
model has antiferromagnetic (AF) long-range order at 
all values of the ratio U/t. For large U/t, each site of the 
lattice is singly occupied, and well defined moments ex- 
ist. AF order arises as a result of a second-order lowering 
of energy when neighboring electron spins are antiparal- 
lel. In this strong-coupling regime, the density of states 
A/'(w) consists of upper and lower Hubbard bands, sepa- 
rated by U. The compressibility k = Nd{n) / vanishes 
at half-filling, reflecting the presence of a Mott-Hubbard 
gap. 

At weak coupling, AF order is produced by nesting 
of the Fermi surface, that is, e(k + Q) = — e(k), at Q = 
(tt, tt), which results in a divergence of the noninteracting 
magnetic susceptibility. Here, e(k) = —2t{coskx -|-cosfcy) 
is the free-particle dispersion relation in the clean limit. 
The density of states exhibits a Slater gap at half-filling, 
arising from this AF order, and again k vanishes. 

Previous DQMC simulations have confirmed this pic- 
ture of the physics of the clean Hubbard model at half- 
filling, and made these statements more quantitative. 
P,0 An analysis of the effect of bond disorder, which we 
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shall review below, has also been performed. [|2j How- 
ever, for the case of site disorder, DQMC simulations 
have not proven possible. In Sec. ||, we review the mean 
field treatment of the problem and consider the effects 
of disorder in this limit. The CPQMC algorithm is out- 
lined in Sec. 



Ill 



^ ^ The effects of disorder on the mag- 
netic correlations are presented in Sec. |^ and we close 
with a brief summary of our results in Sec. 0. The Ap- 
pendix presents some detailed tests of CPQMC on dif- 
ferent model systems (both clean and disordered), with 
a particular focus on the effect of different choices of the 
trial wave functions. 



II. MEAN FIELD APPROXIMATION 

Mean field (MF) theory provides a useful starting point 
for the analysis of the phase diagram of the Hubbard 
model and, as we shall see, also provides us with candi- 
date trial wave functions to use in the CPQMC simula- 
tions. In this approach, the interaction term is decou- 
pled so that the electrons on site i of one spin species 
"see" only the average of the density of the other spin 
species. For zero disorder, the MF phase diagram of the 
2D Hubbard model, as a function of filling and U /t, has 
been given by Hirsch. |13| In the presence of randomness, 
it is important to consider an unrestricted Hartree-Fock 
ansatz (UHF) that allows for general site-dependent oc- 
cupations of each of the two spin species. |Q Systems 
with electron-phonon interactions with quenched lattice 
distortions have been studied in this approximation, 
as have the 3D Anderson-Hubbard model, and the 
propensity for spontaneous phase separation, stripe for- 
mation, and other inhomogeneous charge distributions 
in the clean Hubbard and related models. ||l5| , |l7| , |l8| , ^ 
One of the purposes of this work is to see how such UHF 
results compare to those using CPQMC. 

We will study the disordered 2D model in the UHF 
limit, treating bond and site disorder separately. In order 
to capture the ground state of the model, our calculations 
are performed dX (3 — 1/T — 100, where ks — 1. It is 
useful to define several order parameters for the different 
possible phases. The magnitude of the z component of 
the local moment. 
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measures the tendency for sites to have different numbers 
of up and down spin electrons. The staggered magneti- 
zation. 
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determines the degree of long range antiferromagnetic 
correlation of these moments. Here the notation (. . .) 
represents an averaging over disorder. 



It is also useful to look at charge correlations. Two 
different types of metal-insulator transitions (MIT) can 
occur in the Anderson-Hubbard model. In the Anderson 
MIT, the vanishing of the conductivity is driven primarily 
by the localizing effect of disorder. A useful observable 
is the inverse participation ratio, i?^^. 



R- 
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where the eigenstates of the Hamiltonian read l^ka) = 
■0k(T,i|iCT). For delocalized states, we have Vkcr^ ~ 
\/^/N and Wiajq R^^ ^ in the thermodynamic 
limit. Meanwhile for localized states, the fermions spread 
over a few sites and hence R~^ goes to a finite value in 
the thermodynamic limit, signaling localized electrons. 

In the Mott MIT, the particles are localized primar- 
ily by their interactions. The insulating state is marked 
by the presence of a gap in the charge spectrum at the 
Fermi surface and an associated vanishing of the charge 
compressibility. 



part / 
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where A'part is the total number of particles in the system, 
i.e., iVpart = N{n). 

In related work in 3D, ||l^ Mi has been used to dis- 
tinguish a "spin-glass-like" and a "paramagnetic" disor- 
dered phase that both have Ms = but have Mi nonzero 
and Ml zero, respectively. Interestingly, in two dimen- 
sions, we found the very simple result that Mi = 
whenever Mg = 0. That is, in a MF treatment of the 
half-filled two dimensional model it appears that there 
is no phase in which local moments are present with- 
out ordering antiferromagnetically. However, a spin-glass 
phase has been observed away from half-filling in a sim- 
ilar two dimensional model applicable to the study of 
La2-3;Sr2;Cu04. |2^ As we shall comment further below, 
we are unable to address the possible spin-glass physics 
with QMC, and so we merely report here that, appar- 
ently, within MF theory (MFT) the spin-glass phase is 
absent in two dimensions. Monte Carlo studies of clas- 
sical spin glasses have shown that, as with many phase 
transitions, the appearence of spin-glass order is made 
less likely as the dimensionality is reduced. p^ , |23t 

A similar difference between two and three dimensions 
Jl6t is manifest in the behavior of the inverse partici- 
pation ratio. We have evaluated R~^ but find that it 
approaches finite values (indicating insulating behavior) 
throughout the phase diagram of the disordered two di- 
mensional Hubbard model at half-filling. Since R~^ is 
finite and Mi mimics M^ throughout our phase diagram, 
we will present results only for the staggered magnetiza- 
tion and compressibility. 

Figures. |^ and || show these observables for sweeps of 
the disorder strength at fixed values of the interaction. In 
the case of bond disorder (Fig. ^, the staggered magneti- 
zation is nonzero at all but the smallest interactions. 
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The compressibility, however, has an interesting change 
in behavior as U increases to U = 2t, namely a transition 
from an Anderson insulating phase with nonzero com- 
pressibility to a Mott insulating phase with k — 0. (24j 

For site disorder there is a much clearer region of van- 
ishing staggered magnetization and hence paramagnetic 
behavior. At U = 2t, for example, Ms vanishes beyond 
Vs = At. At U = At, however, the physics becomes rather 
similar to the bond disordered case, with AF correlation 
extending to very large values of disorder, and a signature 
of a Mott transition in the compressibility. The enlarge- 
ment of the paramagnetic phase space at the expense 
of the antiferromagnetically ordered Mott and Anderson 
phases is presumably a consequence of the existence, with 
site disorder, of potential wells on which pairs can form, 
destroying the moments. Similarly, we observe that on- 
site disorder is more effective at eliminating the charge 
gap than bond disorder. 

We used these, and other, sweeps of disorder strength 
at fixed interaction to generate the full UHF ground 
state phase diagrams of the site- and bond-disordered 
two dimensional Hubbard models as shown in Figs. ^ 
and ^. We observed both antiferromagnetically ordered 
Anderson and Mott insulating phases, and a paramag- 
netic Anderson insulating region. These three phases are 
described by 

• Paramagnetic Anderson insulator (P AT): Mi = 0, 

• Antiferromagnetic Anderson insulator (AF AI): 
Mi^O,Ms^O,K^O,R-'^0. 

• Antiferromagnetic Mott insulator (AF MI): Mi ^ 
0, Ms ^0,K = 0,i?-i 0. 

In the case of bond disorder, systems of size 6x6, 
8x8, and 10 x 10 were simulated with averages from 40, 
50, and 50 disorder realizations, respectively. AF order 
dominates the MF phase diagram, as might be expected 
since the disorder is not destroying the moments directly, 
and MFT is too primitive to pick up subtle effects such 
as destruction of AF order via singlet formation. The 
paramagnetic region is restricted to a narrow wedge with 
Vt > 16U. For U >2t = W/4 (where W is the bandwidth 
of the 2D tight binding model) we observe AF ordered 
phases both of the Mott and Anderson variety with a 
boundary given roughly by U = 2t + Vf/4. For U <2t = 
W/A there is no Mott gap within MFT. 



0.75 
Ms 
0.5 

0.25 




8 12 16 



FIG. 1. Staggered magnetization and compressibility for 
the case of bond disorder on a 10 x 10 lattice in the unre- 
stricted Hartree-Fock approximation. AF long-range order 
(LRO) is destroyed by bond disorder for weak on-site inter- 
actions, U ^ 1.0. The staggered magnetization, Ms, levels off 
in the thermodynamic limit for U ~ 2.0, and no amount of 
disorder destroys AF LRO. The inset shows the behavior of 
the compressibility with Vt- A gap is present at U > 2.0 and 
Vt < 2.0, but it is destroyed with increasing disorder. 
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FIG. 2. Evolution of the staggered magnetization and com- 
pressibility for on-site disorder on a 12 x 12 lattice. While 
there is a transition to a paramagnetic state for U = 2, the 
system remains antiferromagnetic even at large 14 for (7 = 4. 
The compressibility (inset) indicates that the presence of a 
gap is more sensitive to site disorder than to bond disorder. 



In the case of site disorder, the antiferromagnetic Mott 
insulator exists for interaction strengths that obey U > 
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Vs- When U > 2t the antiferromagnetic Mott insulator 
is first supplanted by an antiferromagnetic Anderson in- 
sulator with increasing disorder and then ultimately by 
a paramagnetic Anderson insulator at much larger Vg. 
For U < 2t the transition from antiferromagnetic Mott 
insulator appears to go directly to the paramagnetic An- 
derson insulator. 
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FIG. 3. Phase diagram of the bond-disordered Hubbard 
model within the unrestricted Hartree-Fock limit. P = para- 
magnet, AF — antiferromagnet, AI = Anderson insulator, MI 
— Mott insulator. 
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FIG. 4. Phase diagram of the site-disordered Hubbard 
model within the unrestricted Hartree-Fock limit. P — para- 
magnet, AF — antiferromagnet, AI = Anderson insulator, MI 
— Mott insulator. 

In summary, a few notable results from our mean field 
calculation are the following: (i) The same three phases, 
AF ordered Mott insulator, AF ordered Anderson insula- 
tor, and paramagnetic Anderson insulator, are observed 



for bond and site disorder, (ii) We saw no evidence 
for metaUic = 0) or spin-glass-like (Af/ ^ with 

Ms = 0) behavior, (in) AF LRO is never destroyed by 
disorder even at relatively modest values of the interac- 
tion, e.g., J7 ~ 2t for bond disorder and U ^ 4t for site 
disorder. 



III. DESCRIPTION OF THE CPQMC 
SIMULATION 

The ground-state CPQMC method was developed 
to study correlated lattice electrons where no special 
particle-hole symmetry exists to eliminate the sign prob- 
lem. It applies techniques that are a hybrid of determi- 
nant Quantum Monte Carlo (DQMC) and diffusion 
Monte Carlo || (DMC) methods. 

Like the DQMC technique, the method employs the 
Hubbard- Stratonovich (HS) transformation to decouple 
the interaction term of the Hamiltonian, J7 ni|7T.i|. 
The result is a quadratic Hamiltonian in which the in- 
teraction between electrons has been replaced by the in- 
teraction of independent electrons with a classical fluctu- 
ating field. Sampling over the possible values of the HS 
field reproduces the original electron-electron interaction. 

This quadratic Hamiltonian can then be used in an 
imaginary time propagation of a Slater determinant, al- 
lowing projection of the ground state from an initial trial 
wave function: |-0o) — limT-^oo exp(— ri/)|?/'*-°-'). The 
similarity to the DMC method comes from the fact that 
the imaginary time propagation is represented by an en- 
semble of random walkers. (However, the random walk 
in the CPQMC method is performed in a space of Slater 
determinants, in contrast to the DMC method where the 
random walk is in configuration space.) The constrained 
path approximation, necessary for dealing with negative 
weights, is similar in spirit to the fixed node approx- 
imation commonly used to study correlated fermions in 
the continuum. It imposes a boundary condition in de- 
terminant space with a trial wave function, which con- 
strains the random walkers to half of the over-complete 
determinant space. The details of the CPQMC algorithm 
have been discussed elsewhere [Q, but we will provide 
a brief description followed by a discussion of the neces- 
sary adjustments to treat disorder and the observables of 
interest. 



A. Generation of configurations 

At any time in the CPQMC simulation, the wave func- 
tion is represented by an ensemble of random walkers. 
More specifically, we work in a single Slater determinant 
basis and represent our wave function at imaginary time 

step rt by IV'^"^) oc J2k I'^i"'')' where |^^"^) is an individual 
walker (Slater determinant). The initial wave function 
IV'*'"') can in principle be any linear combination of Slater 
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determinants not orthogonal to the ground state. For 
convenience, we choose it to be the unrestricted Hartree- 
Fock wave functions discussed earlier, i.e., IV'*'"') = IV't)- 
In order to propagate the wave function forward to 
imaginary time r, we discretize the propagator to a series 
of short time steps, At. This allows us to apply the 
Trotter approximation 

AtK AtK 
exp(— AriJ) = exp( — ) exp(— AtVF) exp( — ), 

and isolate the potential energy W from the kinetic en- 
ergy K. The HS transformation is then applied. 



exp(— Ar?7ni|nij) — exp 



!;i=±l 



^ii)], (6) 



where cosh(7) = exp(Ar?7/2), p{xi) = 1/2, and Xi is 
the HS field at site i. At each imaginary time step the 
interaction part of the propagator is now a function of 
the HS field. 

Because of the special form of the propagator, each 
Slater determinant in the representation of the 

wave function at imaginary time step n is transformed 
into another Slater determinant 1'/'^"^^'') at imaginary 
time step n -I- 1. We thereby apply the incremental pro- 
jection operator repeatedly to the wave function of our 
system to project out the ground state. In order to make 
the sampling more efficient we can employ an importance 
function to modify the original probability distribution, 
as has been discussed. [|lO| As in the DQMC algorithm, 
the computation time in the CPQMC algorithm scales 
roughly as N^N^^ per random walk step, where N is the 
number of spatial sites in the lattice and iVu, is the num- 
ber of walkers. 

The CPQMC algorithm is exact up to this point. A re- 
maining issue is the constraint to deal with the sign prob- 
lem, which is usually implemented with the importance 
function. We define the overlap integral Ot = {''pT\4'k), 
and demand that individual walkers maintain a positive 
Ot, i.e., that they do not cross the boundary {ipT\4'k) — 
in their random walk in Slater determinant space. This 
is applied to every walker at every time step. The con- 
straint is an approximation, whose quality depends on 
the quality of the trial wave function lipr)- 



B. Measurements 

Monte Carlo methods such as the CPQMC technique 
that employ trial wave functions are well suited for cal- 
culating the ground-state energy, and initial work on the 
CPQMC method demonstrated an excellent agree- 
ment of the ground state energy with exact approaches 
in cases where exact results were available. |^ In the Ap- 
pendix, we will demonstrate that this agreement, though 
a bit less accurate, extends to our simulations. 



In the body of the paper, however, we will focus on real 
space magnetic order which can be identified by measur- 
ing the correlation function. 



i^^(TOjmj+i). 



(7) 



Here toj — — nji is the z component of the lo- 
cal spin operator, and N is the total number of lattice 
sites. C(0, 0) measures the squared local magnetic mo- 
ment (to?). In the clean system, C(0,0) = 0.5 in the 
noninteracting limit, and saturates at C(0, 0) = 1, as C/ 
is increased. In the clean system, (TOjTOj-|_i) is transla- 
tionally invariant, that is, independent of j. For a par- 
ticular disorder realization, however, this is not the case, 
and translational invariance is restored only after disor- 
der averaging. 

It is useful to consider the magnetic structure factor, 
the Fourier transformation of C(l), 



^(q)=^C(l)e 



zql 



(8) 



The structure factor will have sharp peaks at ordering 
vector Q when long-range magnetic order is present. 
(38(0, 0) is the uniform spin susceptibility. At half-filling, 
we always find S'(q) to be largest at the commensurate 
vector Q = (tt^tt), even in the presence of randomness. 
However, our resolution in momentum space is rather 
coarse and ordering at Q values close to (tt, tt) would 
be difficult to see unless the lattice sizes were much in- 
creased. 

For finite lattice simulations, the issue of the presence 
of long-range order in the thermodynamic limit may be 
settled by examining the scaling properties on lattices of 
different size. Spin wave theory pq] predicts 



C{L/2,L/2) = ^ + 0{L-^), 



5'(7r,7r) _ Ml 
N ~ 3 



0{L- 



(9) 



Here Ms is the sublattice magnetization in the thermo- 
dynamic limit, and (L/2,i/2) is the maximal separation 
on a square lattice of linear size L = y/N with periodic 
boundary conditions. C and S provide two quantities to 
extrapolate the value of the ground state order parame- 
ter. 

It is important to comment on the differences of be- 
havior between the order parameters based on local mo- 
ments such as Ml and Ms [Eqs. ||and D, which one might 
expect to find in comparing mean field and QMC treat- 
ments. First, because of the relatively modest spatial 
lattice sizes being simulated in the QMC scheme, over 
the course of a typical run the simulation is able to ex- 
plore the equivalent states that have a surplus of up-spin 
electrons and a surplus of down-spin electrons on a given 
site. In the thermodynamic limit, like a real material, 



5 



the simulation could get stuck in one or the other, but 
our lattice sizes are not large enough for that to happen. 
As a consequence, any direct measurement we make of 
Ms always vanishes, and Ms can only be inferred from 
the finite size scaling analysis, Eq. ^ Meanwhile, MFT 
is able to study Ms directly. Second, MFT has a well- 
known tendency to over-estimate the sharpness of the 
behavior of the local moment as a function of the ratio of 
interaction strength to band-width. For example, even 
in the clean system, phase transitions associated with 
magnetic long-range order are often associated also with 
abrupt formation of local moments. However, it is well 
established that within the QMC method the evolution of 
the local moment is much less abrupt. 1^ For both these 
reasons, we do not expect to be able to observe any spin- 
glass (SG) phase transition within the QMC method, nor 
indeed do we see one. In this work we are only able to 
report that the SG transition observed in MFT at half- 
filling in the three-dimensional model apparently does 
not occur within the same MF treatment in two dimen- 
sions. Spin-glass order in interacting Fermi systems has 
generated much interest recently, but most investigations 
have relied on field theoretic and renormalization-group- 
type techniques. [ p9[ 

The effect of disorder on the size or existence of the 
Mott gap in the CPQMC method is an interesting ques- 
tion that will not be considered in detail in this work. 
We have measured the density as a function of chemical 
potential |Q and the compressibility, and find the Mott 
gap is rather strongly reduced by site disorder and con- 
siderably less aff'ected by bond disorder, but we leave 
a detailed analysis to a later presentation. 



IV. RESULTS FOR MAGNETIC CORRELATIONS 

In this section we address the primary point of interest 
in the paper, the effect disorder has on the long-range 
magnetic order of the half-filled 2D Hubbard model. 
In particular, we want to determine the critical disor- 
der strength necessary to destroy the magnetically or- 
dered ground state and ascertain the accuracy of the 
UHF phase diagram. We concentrate on the case where 
U — At, where the UHF shows no transition to a param- 
agnetic order. This antiferromagnetic to paramagnetic 
transition can only take place in the thermodynamic 
limit; hence we resort to finite size scaling techniques 
to calculate the disorder strength at which AF LRO is 
lost. In the case of bond disorder, these questions have 
been investigated previously with the DQMC method 
since there is no sign problem. We re-address these 
questions here in order to benchmark the accuracy of the 
CPQMC method. After the results for bond disorder are 
discussed, we consider the case of random site energies 
and its ability to drive the model from an antiferromagnet 
to a paramagnet. We again point out that the site dis- 
ordered problem has not yet been studied with th QMC 



method because of the sign problem. 



A. Bond disorder 

The mechanism by which bond disorder destroys AF 
LRO is the formation of spin-0 singlets. When the dis- 
order becomes large, the lattice contains strong bonds 
where it is favorable for nearest neighbors to form sin- 
glets rather than participate in AF LRO. Unlike the case 
of site disorder, one would expect the persistence of local 
correlations, C(0,0), even in the paramagnetic phase. 

In our simulations of this problem, we draw random 
near-neighbor hopping strengths from a uniform dis- 
tribution about a mean value of (iy) = t — 1, e.g., 
[1 — Vt/2,1 -|- Vt/2]. We considered simulations with a 
renormalized interaction, C/cpqmc U , chosen, as de- 
scribed in the Appendix, so the CPQMC and DQMC 
results match in the clean limit. For 4 x 4, 6 x 6, and 
8x8 lattices the values of ?7cpqmc were set to 1.25, 1.75, 
and 1.85, respectively. Disorder averaging was performed 
on 10 realizations for each disorder strength on each lat- 
tice. We used as \'4't) an UHF trial wave function (TWF) 
obtained with ?7twf = 2, but as shown in the Appendix, 
the results are not sensitive to this choice. |Q 

The real space spin correlations are shown in Fig. || on 
a fixed 8x8 lattice size. Disorder substantially decreases 
the antiferromagnetic order, with only a relatively minor 
suppression of the squared local moment (inset). 

Summing up these real space spin correlations yields 
the magnetic structure factor. Working on a range of 
lattice sizes, and employing the scaling relationship for 
5'(7r, tt), Eq. yields the staggered magnetization as a 
function of disorder from the intercepts of our scaling 
plot. Fig. I 




FIG. 5. The real space spin-spin correlations on a 
bond-disordered 8x8 lattice. The inset shows the squared 
local moment scaled by the noninteracting value as a function 
of bond disorder strength. 
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FIG. 6. Finite size scaling for the structure factor as 
a function of bond disorder. Simulations were performed 
with a renormalized interaction and an UHF trial state with 
f/TWF = 2. Extrapolations to the thermodynamic limit 
give an intercept that is equal to /Z. A critical disor- 
der strength of Vt ~ 1.6 ± 0.4 is found, which agrees with the 
DQMC result. The dashed lines are linear least squares fits 
to the data. 



ture factor S'(7r, tt) were done. |^ The error bars on C 
were only slightly larger than those from 5, so they pro- 
vided essentially equivalent (and consistent) information. 
In the disordered system, however, we have found that 
C is not as useful to examine as the structure factor. We 
believe that the reason is that on a lattice of N sites, 
C is obtained through an average over only the 0{N) 
pairs of sites separated by (L/2, L/2) while S involves all 
0{N'^) separations. The structure factor is thus much 
less sensitive to individual disorder realizations. Second, 
the linear extrapolation to negative order parameter in 
the disordered phase evident in Fig. || has been observed 
in previous simulations, e.g., of the clean periodic An- 
derson model as the hybridization between localized and 
delocalized orbitals is increased until the antifcrromag- 
netism gives way to Kondo singlet formation, psf The 
point is that Eq. ^, which yields a linear extrapolation in 
1/L, is valid only in the ordered phase. In the disordered 
phase S/N cx b/L^ so one should actually fit the data to 
a quadratic expression that goes through the origin. 

The results for Ms are exhibited in Fig. 0, which gives 
a critical disorder of Vt = 1.6 ±0.4. The UHF calculation 
predicts no transition to a paramagnetic phase for this 
value of U, but the CPQMC calculation agrees very well 
with the previous DQMC results. We did not observe 
an enhancement of at weak disorder as was observed 
in our UHF data and in DMFT and DQMC calculations, 
1^ , ^ but our resolution might be too small to observe 
this effect. 
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FIG. 7. Staggered magnetization as a function of bond dis- 
order. Values were calculated from the intercepts of Fig. ^ 
Ms vanishes at the critical disorder strength Vt « 1.6 ± 0.4. 

It is worth making two further comments conncerning 
Fig. 6. First, in papers establishing long-range antifcrro- 
magnetic order in the clean, half-filled, two-dimensional 
Hubbard model, extrapolations of both the longest-range 
spin-spin correlation function C{L/2, L/2) and the struc- 



B. Site disorder 

In our simulations of the site-disordered model, ran- 
dom energies were selected from a uniform distribution 
Ei e [-~Vs/2,Vs/2]. Sites with ei < favor double occu- 
pancy while sites with ej > favor the unoccupied state. 
This leads to a direct destruction of moments, unlike the 
case of bond disorder. In the presence of a repulsive Hub- 
bard interaction U , there is therefore a competition be- 
tween a lattice with local moments and AF LRO, which is 
favored by ?7, and a state of doubly occupied and empty 
sites favored by the disorder. 

As discussed in the Appendix, simulations with on-site 
disorder need to have C/twf > Vs in order to capture 
the physics of the model and not the effect of trial wave 
function. We used a trial wave function with J7twf = 6, 
and the same renormalized interaction J/cpqmc used in 
the bond-disordered case. We simulated 10 reahzations 
of each disorder strength. The suppression of the real 
space spin correlations is displayed in Fig. ^. 

We can again analyze appropriately scaled data on dif- 
ferent lattices, with the results shown in Fig. |^. The 
intercepts give the staggered magnetization, which is 
driven to zero above a critical site disorder strength of 
5.0 ± 0.5 (Fig. 



no transition to a paramagnetic phase for this value of U . 



10). The UHF calculation predicts 
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Hence, our UHF trial wave function has long-range an- 
tiferromagnetic correlations throughout the range of Vs 
in Fig. |l^, and the destruction of order is not associated 
with any change in the nature of \iPt)- 
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FIG. 8. The effect of site disorder on the real space 
spin-spin correlations on an 8 x 8 lattice. Correlations for 
distances greater than 1 are uniformly reduced by disorder. 
The inset shows the behavior of the scaled squared local mo- 
ment as a function of Vs. Site disorder strongly suppresses 
the local moment. 



0.25 



0.2 



.0.15 



W 0.1- 



0.05 





=0.0 


V: 

S 


=1.0 


v = 

s 


=2.0 




=3.0 


V: 

S 


=4.0 




=4.5 


v = 

s 


=5.0 


V: 

S 


=5.5 




=6.0 



0.05 



0.1 



0.15 



0.25 



0.3 



1/L 



FIG. 9. Scaling relationship for the site disordered model 
with a renormalized (7cpqmc and an UHF trial state with 
Utwf ~ 6. The critical disorder strength was Vs ~ 5.0 ± 0.5. 
The dashed lines are linear least-squares fits to the data. 




V. 



FIG. 10. Staggered magnetization (A/s) for the 
site-disordered 2D Hubbard model. Data were obtained from 
the intercepts of the scaling relationship for S{n,n), Fig. ^. 
Ms vanishes at Vs ^ 5.0 ± 0.5. 



0.45 

0.4 
0.35 

0.3 
Q0.25 

0.2 
0.15 

0.1 



o 8x8, U =1.85 

' cpqmc 



V. 



FIG. 11. The double occupancy as a function of site dis- 
order on an 8 X 8 lattice. The dashed line denotes the sep- 
aration between the repulsive (below) and attractive (above) 
Hubbard model ai t — and any temperature. Data were ob- 
tained from simulations with ?7cpqmc ~ 1.85 and Utwf = 6. 

It is interesting that the point at which AF LRO is 
lost corresponds rather closely to the value of randomness 
where the squared local moment is reduced below its non- 
interacting value. This is emphasized in Fig. |ll|, which 
shows the average double occupancy D — (n^-ni) = 
Ei/UcpQMC in the CPQMC method, where Ej is the in- 
teraction energy of the fermions. Since C(0, 0) = 1 — 2D, 
an enhancement of D above the noninteracting value is 
synonymous with the moment falling below the U = 
value. 
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V. CONCLUSIONS 

We have studied the ground state half-filled 2D Hub- 
bard model with both bond and site disorder in the mean 
field limit and by the constrained path quantum Monte 
Carlo method. Our most significant quantitative result 
was the first computation of the critical disorder strength 
for the destruction of antiferromagnetic long-range order 
by site randomness, {Vs)ci:it — 5i for U = At. For this 
value of the interaction strength, no amount of site dis- 
order destroys the order in the Hartrce-Fock approach, 
so this emphasizes the need for better treatment of cor- 
relations that techniques like the quantum Monte Carlo 
method provide. In general, we find that UHF calcu- 
lations grossly overestimate the tendency for magnetic 
order when compared to the CPQMC calculations, as 
might be expected by an approach which ignores fluctu- 
ations. 

There is less significant disagreement between UHF 
and CPQMC techniques for the transport properties. 
Unlike the 3D case, the UHF treatment finds that the 
2D Hubbard model is insulating for all values of inter- 
action and disorder at half-filling: the inverse participa- 
tion ratio was always nonzero. This is consistent with 
recent QMC calculations in two dimensions. |3^] A fur- 
ther difference between the 2D and 3D UHF results is 
our conclusion that the spin-glass phase is absent in two 
dimensions at half-filling. It is interesting to note that 
there have been some indications in the QMC treatment 
of a metal-insulator transition off half-filling in the 2D 
Hubbard model. |^ 

Our work further evaluated and extended the range of 
validity of the CPQMC method by applying it to random 
systems. We found that, as is the case for clean systems, 
the CPQMC technique can provide an accurate way of 
treating the Hubbard model. In particular, it gives the 
same critical disorder strength as the DQMC method in 
the case when DQMC has no sign problem (bond ran- 
domness), which provides some confidence in applying 
it to cases such as site-disordered problems where the 
DQMC method cannot give reliable results. 

Much of the initial theoretical evidence for, and under- 
standing of, questions of charge ordering in Hubbard-like 
models has come from UHF treatments, jl^] Recent work 
with techniques such as the density matrix renormaliza- 
tion group have emphasized that stripe formation is a 
subtle and delicate effect. Our work indicates that 
there are significant corrections to the spin correlations 
within UHF treatments, and that further CPQMC calcu- 
lations hold promise to shed some light on the behavior 
of disordered and interacting electron systems. 
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APPENDIX: TESTS OF THE ALGORITHM 

In this Appendix we describe some tests of the 
CPQMC algorithm with a specific focus on the effect of 
the trial wave functions on the results. Previously, the 
CPQMC algorithm has been extensively tested for in- 
teracting systems with no disorder. Investigations have 
been performed on the single-band Hubbard model in re- 
gions of parameter space where a severe sign problem is 
known to exist. The method has also been used to 
the study superconductivity in the 2D Hubbard model by 
looking for long-range pairing correlations in the ground 
state and for ferromagnetism in the periodic Ander- 
son model. p7|] Here disorder has been considered within 
a CPQMC calculation. 

CPQMC is an exact algorithm, for all observables, in 
the absence of interactions, even when disorder is turned 
on. As a check of our code, we therefore first verified that 
the CPQMC code reproduced results from exact diago- 
nalization (ED) at different disorder strengths. Likewise, 
the DQMC algorithm agrees perfectly with ED results. 

We next looked in detail at the behavior of the mag- 
netic structure factor as a function of disorder and in- 
teraction strengths, and as a function of the trial wave 
function. In the following discussion it is useful to dis- 
tinguish between the value of the interaction strength, 
C^cpqmCj used in the CPQMC algorithm, the value of 
the interaction strength, J7twf, used in the trial wave 
function, and the physical value of U in the Hamilto- 
nian. We concentrate here on the case where U = At and 
t = 1. In our DQMC simulations, which we used to pro- 
vide benchmarks for the CPQMC method, we of course 
always choose C/dqmc = U since the DQMC treatment 
is exact. 

The key result of our studies is that the CPQMC tech- 
nique with UHF trial wave functions significantly overes- 
timates the magnetic correlations if ?7cpqmc = U . This 
is illustrated in Fig. which shows the ratio of the 
CPQMC to DQMC structure factors for the clean sys- 
tem on a 4 X 4 lattice as a function of C/cpqmc/C^- The 
results are independent of the value of t/xwF as long as 
t^TWF 7^ 0. The structure factor is the same for the two 
methods when C/cpqmc/C^ ~ 0.31 or C/cpqmc ~ 1-25. 
For 6x6 and 8x8 lattices agreement is attained at 
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t^CPQMC 

[7 = 4. 



1.75 and L'cpqmc ~ 1.85, respectively, for wave function. 



4x4 lattice, half-filled, U.^=4 V^=0.01 

aqmc s 




^cpqmc/^dqmc 
FIG. 12. Scaled data for S(Tr,n)cpQMC / S{-k,tt)dqmc as 
a function of J/cpQMc/t^DQMC for the clean Hubbard model 
in the CPQMC treatment. When (7cpqmc ~ f/DQMc(= U), 
AF LRO is considerably overestimated. We found agreement 
between the two methods at (7cpqmc ~ 1-25 for this 4x4 
lattice. The interaction strength Utwf of the UHF trial wave 
function does not affect results as long as Utwf 0. 

A similar effect is seen when disorder is turned on, as 
illustrated in Fig. |l^ for site disorder. Here we show 
the ratio of the CPQMC and DQMC structure factors 
as a function of f/cPQMc/C^ for different lattice sizes. 
These results are for a single realization of disorder, 
which is kept fixed as C/cpqmc and C/twf are varied. 
The values of the CPQMC structure factor for differ- 
ent Utwf fall onto two curves: For all C/twf > the 
CPQMC method gives the same significantly overesti- 
mated structure factor. Meanwhile, for all C/twf < 
the CPQMC methods gives the same significantly under- 
estimated structure factor. 

These different behaviors are a direct manifestation of 
the trial wave functions. In the unrestricted Hartree- 
Fock calculation, both C/twf and Vs act as one-body 
potentials and compete with each other: When C/twf < 
Vs, Vs dominates and double occupancy is allowed, i.e., 
the system is more free-electron like; when C/twf > Vs, 
C/twf dominates, double occupancy is discouraged, and 
the system prefers to be in an AF state. Clearly, the 
CPQMC technique could not adequately eliminate the 
biases that the two different classes of UHF trial wave 
functions introduce through the approximate constraint 
to bring quantitative agreement between the two sets of 
results. In all the work reported in the body of this paper 
we chose C/twf > Vs so that the trial wave function had 
long-range antiferromagnetic order. The destruction of 
order as randomness increased therefore must occur from 
correlation effects and not from any transition in the trial 
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FIG. 13. Scaled data for S{tv,-k)cpqmc/S{it,-k)dqmc for 
the site-disordered Hubbard model in the CPQMC technique: 
(a) 4 X 4, K = 2; (b) 4 x 4, = 4; (c) 6 x 6, K = 2. 
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Bond disorder is studied in Fig. |T|. The CPQMC 
structure factor is again overestimated. Bond disorder 
Vt, however, does not turn into one-body potentials in 
the UHF method in a simple manner, and unlike site dis- 
order, there are not two separate behaviors. The result 
is rather independent of C/twf- 
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FIG. 14. Scaled data for S{-k,tt)cpqmc / S{-!t,tt)dqmc for 
the bond-disordered Hubbard model in the CPQMC method: 
(a) 4 X 4, K = 2; (b) 6 x 6, 14 = 2. 

While we show results for single disorder realizations in 
Figs. |l3| and |l^, the values of the renormalized couplings 
used in determining, for example, the critical disorder 
strength, in the main text of this paper were obtained 
with disorder averaging. The overestimation of the struc- 
ture factor is a significant concern in our studies of the 
destruction of long-range antiferromagnetic order, which 
rely primarily on this quantity. It exemplifies the diffi- 
culty that faces all approximate methods to deal with the 
sign problem that use a trial wave function to constrain 
the QMC sampling, namely the results can be biased by 
the trial wave function, sometimes significantly. Our ap- 



proach is to fix C/cPQMC at a "renormalized" value so 
that the structure factor from the CPQMC algorithm 
matches the DQMC result. This sort of tuning of the 
interaction strength has previously been done in com- 
parisons of diagrammatic calculations for the Hubbard 
model with DQMC results. |^ A crucial question, of 
course, is whether the renormalized J7cpqmc is indepen- 
dent of lattice size. We found that J/cpqmc depends only 
weakly on lattice size for L > 4, as seen in Fig. 1^. Again, 
similar effects are known in the comparisons of DQMC 
and diagrammatic calculations. pg| ] 

A further indication of the importance of the renor- 
malization of the interaction lies in the behavior of the 
staggered magnetization per site. Data using a renormal- 
ized ?7cPQMC always lie below the classical upper limit 
of 0.5 whereas simulations for fixed C/cpqmc = U did 
not. At Vt — 0, our result with a renormalized ?7cpqmc 
and C/twf = 2 was Alg — 0.33(2). This clean sys- 
tem value compares well to earlier results for the quan- 
tum Heisenberg model obtained from a QMC calcula- 
tion, 0.30(2), |40| and from perturbation series expan- 
sions, 0.313. |§ 

Another crucial question is whether the renormalized 
C^CPQMC depends on disorder strength. This question can 
be addressed in the case of bond disorder where DQMC 
simulations of large lattices at low T can be done with- 
out encountering the sign problem, but cannot be done 
for site disorder. We found that for a given lattice size 
a single constant choice of C/cpqmc could be used for 
all Vt- We note that the apparent variation of the 
renormalization on the values of Vs and Vt evident in 
comparing Figs. 12-14 is dominantly due to the fact the 
data presented there are not disorder averaged, a particu- 
larly important issue for the smaller 4x4 lattices. When 
such averaging is done, as in the main body of this paper, 
the variation is very significantly reduced. This is fortu- 
nate, since the tuning of C/cpqmc for different disorder 
strength would be not only awkward but would also call 
into question whether transitions we observe as a func- 
tion of disorder strength were caused by the tuning or by 
the randomness itself. 

We have focused here on the behavior of the struc- 
ture factor and matching the DQMC and CPQMC val- 
ues. Previous work has shown that the energy and other 
correlations agree well. [ p^|j3^j37| | We have verified that 
the energy in DQMC and CPQMC techniques remains 
in relatively good agreement in these disordered systems 
if the trivial difference in interaction energy is accounted 
for by defining, 



rprn 

^CPQMC 



EcPQMC + {U-DQMG — UcPQMc){D 



(10) 



Comparisons of the energy in CPQMC and DQMC tech- 
niques behave as shown in Fig. na. For the parameters 
used in our simulations, the energies disagree by at most 
5%. 
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